home *** CD-ROM | disk | FTP | other *** search
/ Technotools / Technotools (Chestnut CD-ROM)(1993).ISO / lang_c / cug167 / spline.doc < prev   
Text File  |  1985-08-21  |  11KB  |  304 lines
  1. /* 
  2. HEADER:     CUG
  3. TITLE:        Cubic Spline Functions Theory
  4. VERSION:    1.00
  5. DATE:        11/01/85
  6. DESCRIPTION:    "Mathematical background and development of
  7.         equations used in SPLINE, a full emulation of
  8.         Unix's 'spline' utility."
  9. KEYWORDS:    spline, graphics, plot, filter, UNIX
  10. SYSTEM:
  11. FILENAME:    SPLINE.DOC
  12. WARNINGS:
  13. CRC:        xxxx
  14. SEE-ALSO:    SPLINE.C
  15. AUTHORS:    Ian Ashdown - byHeart Software
  16. COMPILERS:
  17. REFERENCES:    AUTHORS: R.W. Hamming
  18.         TITLE:     Numerical Methods for Scientists and
  19.              Engineers, 2nd Edition
  20.              pp. 349 - 355, McGraw-Hill (1973)
  21.         AUTHORS: Forsythe, Malcolm & Moler
  22.         TITLE:     Computer Methods for Mathematical
  23.              Computations
  24.              pp. 70 - 83, Prentice-Hall
  25. ENDREF
  26. */
  27.  
  28. /*-------------------------------------------------------------*/
  29.  
  30. byHeart Software
  31. 1089 West 21st Street
  32. North Vancouver, B.C. V6J 1G8
  33. Canada
  34. November 1st, 1985
  35.  
  36.           Cubic Spline Functions Theory
  37.           -----------------------------
  38.  
  39.              Ian E. Ashdown
  40.             byHeart Software
  41.  
  42. References:    1. Numerical Methods for Scientists and Engineers
  43.            R.W. Hamming, 2nd Edition
  44.            pp.349 - 355, McGraw-Hill (1973)
  45.         2. Computer Methods for Mathematical Computations
  46.            Forsythe, Malcolm & Moler
  47.            pp. 70 - 83, Prentice-Hall
  48.  
  49.     Cubic spline functions are a mathematical abstraction of a
  50. draftsman's mechanical splines. Given a set of abscissae and
  51. ordinates, the draftsman first plots this data on a sheet of
  52. paper, then uses a spline - a flexible strip of wood or plastic -
  53. to draw a smooth curve between them. The spline is held in place
  54. at the given points (historically called "knots") with weights.
  55. If the spline is not too severely stressed, elementary beam
  56. theory shows that it will conform to a curve described by a cubic
  57. polynomial equation between each pair of knots, and that adjacent
  58. polynomials will join continuously with continuous first and
  59. second derivatives. The result is a visually "smooth" curve. The
  60. following discussion develops the mathematics necessary to model
  61. these splines with a computer.
  62.  
  63.     Given a set of datum, stated as abscissae x[i] (i = 1,...,n)
  64. and corresponding ordinates y[i], let y(x) be an interpolated
  65. curve through this data, defined as the composite function
  66.  
  67.     y(x) = f[i](x) for x[i] <= x < x[i+1], i = 1,...,n-2  (1)
  68.                and x[n-1] <= x <= x[n]
  69.  
  70. where each function f[i](x) is a cubic polynomial of the form
  71.  
  72.     f[i](x) = a * pow(x,3) + b * pow(x,2) + c * x + d     (2)
  73.  
  74.     (a,b,c and d are constants, and "pow(x,y)" means "x" to
  75.     the power of "y")
  76.  
  77.     In other words, f[](x) is a set of functions, each of which
  78. is defined only over the specific interval between two adjacent
  79. abscissae. (This is equivalent to describing the spline between
  80. each pair of adjacent knots using a cubic polynomial equation, as
  81. per elementary beam theory.) The concatenation of these functions
  82. over the given range of abscissae defines the function y(x).
  83.  
  84.     Furthermore, let's define y'[i] and y"[i] as the first and
  85. second derivatives of y(x) at x = x[i]. Continuous first and
  86. second derivatives at the knots of the spline then imply the
  87. following continuity conditions:
  88.  
  89.     f[i](x[i]) = y[i]        i = 1,...,n-1          (3)
  90.     f[i-1](x[i]) = y[i]        i = 2,...,n          (4)
  91.     f'[i-1](x[i]) = f'[i](x[i])    i = 2,...,n-1          (5)
  92.     f"[i-1](x[i]) = f"[i](x[i])    i = 2,...,n-1          (6)
  93.  
  94.     These conditions state that the functions f[i](x) are
  95. continuous where they join at their endpoints (the ordinates
  96. y[i]), thus ensuring the "smoothness" of the composite curve
  97. y(x).
  98.  
  99.     Since each function f[i](x) is a cubic polynomial, it follows
  100. that its third derivative is a constant. This means that y"(x) is
  101. a linear function over each interval. Thus, defining
  102.  
  103.     h[i] = x[i+1] - x[i]                      (7)
  104.  
  105. linear interpolation gives the equation
  106.  
  107.     f"[i](x) = (y"[i] * (x[i-1] - x) +               (8)
  108.            y"[i+1] * (x - x[i]))/h[i]
  109.  
  110.     Integrating this equation twice and selecting the constants
  111. of integration such that Equations (3) and (4) are satisfied
  112. gives:
  113.  
  114.     f[i](x) = (y[i] * (x[i+1] - x) +              (9)
  115.           y[i+1] * (x - x[i]))/h[i] -
  116.           (pow(h[i],2)/6 *
  117.           (y"[i] * ((x[i+1] - x)/
  118.           h[i] -(pow((x[i+1] - x)/h[i],3)) -
  119.           y"[i+1] * ((x - x[i])/h[i] -
  120.           pow((x - x[i])/h[i],3)))
  121.  
  122.     Note that the first two terms are a linear interpolation
  123. function between ordinates y[i] and y[i+1]. The remainder of the
  124. equation is the cubic "correction factor".
  125.  
  126.     Differentiating Equation (9) and evaluating yields:
  127.  
  128.      f'[i](x[i]) = (y[i+1] - y[i])/h[i] - h[i] *         (10)
  129.               (2 * y"[i] + y"[i+1])/6
  130.  
  131.      f'[i-1](x[i]) = (y[i] - y[i-1])/h[i-1] + h[i-1] *    (11)
  132.             (y"[i-1] + 2 * y"[i])/6
  133.  
  134.     Using the identity shown in Equation (5), we get:
  135.  
  136.     h[i-1] * y"[i-1] + 2 * (h[i-1] +             (12)
  137.     h[i]) * y"[i] + h[i] * y"[i+1] =
  138.     6 * ((y[i+1] - y[i])/h[i] - (y[i] - y[i-1])/h[i-1])
  139.  
  140.     for i = 2,...,n-1
  141.  
  142. which must be satisfied at n-2 points by the n unknown quantities
  143. y"[i]. In matrix form (using n = 6 as an example), we have:
  144.  
  145. +-                   -++-     -+       +-    -+    (13)
  146. | h[1] c[2] h[2] 0    0    0    || y"[1] | = 6 * | d[2] |
  147. | 0    h[2] c[3] h[3] 0    0    || y"[2] |       | d[3] |
  148. | 0    0    h[3] c[4] h[4] 0    || y"[3] |       | d[4] |
  149. | 0    0    0    h[4] c[5] h[5] || y"[4] |       | d[5] |
  150. +-                   -+| y"[5] |       +-    -+
  151.                  | y"[6] |
  152.                  +-     -+
  153.  
  154. where c[i] = 2 * (h[i-1] + h[i])
  155. and   d[i] = (y[i+1] - y[i])/h[i] - (y[i] - y[i-1])/h[i-1]
  156.  
  157.     Note that if we limit the abscissae to being equally spaced,
  158. Equation (12) becomes:
  159.  
  160.     h * y"[i-1] + 4 * h * y"[i] + h * y"[i] =         (14)
  161.     6 * (y[i+1] - 2 * y[i] - y[i-1])/h
  162.  
  163.     for i = 2,...,n-1
  164.  
  165. and our example Equation (13) becomes:
  166.  
  167.     +-         -++-     -+       +-    -+             (15)
  168.     | 1 4 1 0 0 0 || y"[1] | = 6 * | d[2] |
  169.     | 0 1 4 1 0 0 || y"[2] |       | d[3] |
  170.     | 0 0 1 4 1 0 || y"[3] |       | d[4] |
  171.     | 0 0 0 1 4 1 || y"[4] |       | d[5] |
  172.     +-         -+| y"[5] |       +-    -+
  173.                | y"[6] |
  174.                +-     -+
  175.  
  176.     where d[i] = (y[i+1] - 2 * y[i] - y[i-1])/h
  177.  
  178.     Two more conditions are required to obtain a unique solution.
  179. These can be obtained by specifying one condition at each end of
  180. the interpolated curve. Several variations are possible, with
  181. two being examined here.
  182.  
  183.     The first one is to specify that
  184.  
  185.     y"[1] = k * y"[2] and y"[n-1] = k * y"[n]         (16)
  186.  
  187. where k is a constant. Equation (13) then becomes:
  188.  
  189. +-                   -++-     -+       +-    -+    (17)
  190. | -1   k    0    0    0    0    || y"[0] | = 6 * | 0    |
  191. | h[1] c[2] h[2] 0    0    0    || y"[1] |     | d[2] |
  192. | 0    h[2] c[3] h[3] 0    0    || y"[2] |       | d[3] |
  193. | 0    0    h[3] c[4] h[4] 0    || y"[3] |       | d[4] |
  194. | 0    0    0    h[4] c[5] h[5] || y"[4] |       | d[5] |
  195. | 0    0    0    0    k    -1    || y"[6] |     | 0    |
  196. +-                   -++-     -+       +-    -+
  197.  
  198. where c[i] = 2 * (h[i-1] + h[i])
  199.       d[i] = (y[i+1] - y[i])/h[i] - (y[i] - y[i-1])/h[i-1]
  200.  
  201.     For the case of equidistant abscissae, this becomes (with a
  202. bit of matrix manipulation):
  203.  
  204.     +-         -++-     -+       +-    -+         (18)
  205.     | 4-k 1   0   0   || y"[2] | = 6 * | d[2] |
  206.     | 1   4   1   0   || y"[3] |       | d[3] |
  207.     | 0   1   4   1   || y"[4] |       | d[4] |
  208.     | 0   0   1   4-k || y"[5] |       | d[5] |
  209.     +-         -++-     -+       +-    -+
  210.  
  211.     where d[i] = (y[i+1] - 2 * y[i] - y[i-1])/(h * h)
  212.  
  213.     If the value of k is zero, we have y"[1] = y"[n] = 0. This is
  214. equivalent to a mechanical spline whose ends are not constrained
  215. beyond the data endpoints, and is known as the "natural" cubic
  216. spline. A nonzero value of k is equivalent to bending the ends of
  217. the mechanical spline, and will affect all of the interior cubic
  218. polynomial functions. In this sense, the cubic spline is a global
  219. function; however, the effect of changing k on the interior
  220. polynomials rapidly decreases as one moves away from the
  221. endpoints (which is what one would expect from intuition).
  222.  
  223.     The second variation is more interesting, and comes from the
  224. need to interpolate data extracted from periodic phenomenon. A
  225. good example of this is a closed curve such as an ellipse. If the
  226. abscissae are expressed over 360 degrees, the curve can be
  227. modeled by a cubic spline function. However, we need to somehow
  228. specify that the endpoints of the curve must be continuous with
  229. respect to each other.
  230.  
  231.     The ordinates are constrained by the data given to be equal
  232. (i.e. - y[1] = y[n]). We need to specify the end conditions such
  233. that 
  234.  
  235.     y'[1] = y'[n] and y"[1] = y"[n]                 (19)
  236.  
  237. The second derivatives are easy - they can be expressed directly
  238. in matrix form. However, to solve for the first derivative of
  239. y(x) at the end points, we need an equation that relates it to
  240. y(x) and its second derivative. Equations (10) and (11) satisfy
  241. this requirement. Evaluating them for x[1] and x[n] respectively
  242. produces:
  243.  
  244.     f'[1](x[1]) = (y[2] - y[1])/h[1] - h[1] *         (20)
  245.               (2 * y"[1] + y"[2])/6
  246.  
  247.     f'[n-1](x[n]) = (y[n] - y[n-1])/h[n-1] + h[n-1] *    (21)
  248.             (y"[n-1] + 2 * y"[n])/6
  249.  
  250. But y'[1] = y'[n], so that
  251.  
  252.     6 * (y[2] - y[1])/h[1] - (y[n] - y[n-1])/h[n-1] =    (22)
  253.     h[n-1] * (y"[n-1] + 2 * y"[n]) +
  254.     h[1] * (2 * y"[1] + y"[2])
  255.  
  256.     Our example equation (13) thus becomes:
  257.  
  258. +-                   -++-     -+       +-    -+    (23)
  259. | 1    0    0    0    0    -1   || y"[0] | = 6 * | 0    |
  260. | h[1] c[2] h[2] 0    0    0    || y"[1] |     | d[2] |
  261. | 0    h[2] c[3] h[3] 0    0    || y"[2] |       | d[3] |
  262. | 0    0    h[3] c[4] h[4] 0    || y"[3] |       | d[4] |
  263. | 0    0    0    h[4] c[5] h[5] || y"[4] |       | d[5] |
  264. | 2a   a    0    0    b    2b    || y"[6] |     | e    |
  265. +-                   -++-     -+       +-    -+
  266.  
  267. where  c[i] = 2 * (h[i-1] + h[i])
  268.        d[i] = (y[i+1] - y[i])/h[i] - (y[i] - y[i-1])/h[i-1]
  269.        a = h[1]
  270.        b = h[n-1]
  271.        e = (y[2] - y[1])/h[1] - (y[n] - y[n-1])/h[n-1]
  272.  
  273.     For equidistant abscissae, this becomes (again with some
  274. matrix manipulation):
  275.  
  276.     +-        -++-     -+       +-    -+         (24)
  277.     | 4 1 0 0 1  || y"[2] | = 6 * | d[2] |
  278.     | 1 4 1 0 0  || y"[3] |       | d[3] |
  279.     | 0 1 4 1 0  || y"[4] |       | d[4] |
  280.     | 0 0 1 4 1  || y"[5] |       | d[5] |
  281.     | 1 0 0 1 4  || y"[6] |       | e    |
  282.     +-        -++-     -+       +-    -+
  283.  
  284.     where d[i] = (y[i+1] - 2 * y[i] - y[i-1])/(h * h)
  285.           e = (y[2] - y[1] - y[n] + y[n-1])/(h * h)
  286.  
  287.     Other end conditions are possible. Using Equations (10) and
  288. (11), you can specify the slope of the splines outside of the
  289. range of the data by specifying the first derivatives at y[1] and
  290. y[n]. You can also specify a linear combination of first and
  291. second derivatives at the end points. However, the above will
  292. generally prove the most useful.
  293.  
  294.     In closing, it's worth noting that the matrices presented
  295. above are diagonally dominant band matrices. This means that the
  296. equations they represent will always have a solution, and that
  297. they can be solved by reduction to upper triangular form and
  298. subsequent back substitution without the need for pivoting. A few
  299. moments' thought will also show that the non-zero elements can be
  300. stored in a few linear arrays, and the remainder ignored. This
  301. means that simple dedicated algorithms can be written to solve
  302. the matrices in linear rather than cubic time, and to do so with
  303. a very efficient utilization of core memory.
  304. 2] -